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Abstract 

We derive an exact recursion formula for the calculation of thermodynamic 
functions of finite systems obeying Bose-Einstein statistics. The formula is 
applicable for canonical systems where the particles can be treated as non- 
interacting in some approximation, e.g. like Bose-Einstein condensates in 
magnetic traps. The numerical effort of our computation scheme grows only 
linear with the number of particles. 

As an example we calculate the relative ground state fluctuations and spe- 
cific heats for ideal Bose gases with a finite numbers of particles enclosed in 

containers of different shapes. 
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With the observation of Bose-Einstein condensation (BEC) of magnetically and 
optically trapped atoms new insights into the nature of this state of matter have been 
given. The experimental situation is in all cases quite different from the ideal gas treated 
within the grand canonical ensemble, which is the standard textbook example. First, the 
number of particles within the traps is fixed and finite, which suggests a canonical or micro- 
canonical treatment of the systems. Second, the confining trap potentials greatly influence 
the condensate properties. Third, although the trapped gases are quite dilute the validity 
of the treatment as non-interacting particle gases has to be checked from case to case. 
Even within the approximation of non-interacting particles the calculation of the thermo- 
dynamic properties of the Bose-Einstein systems remains a hard to tackle mathematical 
problem. Recently, some approximate methods to calculate the fluctuation of the ground 
state occupation number in a trapped Bose-Einstein condensate have been developed ||5|-§ . 
Here we present an exact method to calculate all thermodynamic quantities of finite canon- 
ical Bose systems, given the one particle density of states. 

As the starting point we utilize the recursive formula of the canonical partition function 
for a system of noninteracting bosons as given in 

1 ^ 

ZN{P) = j,Y.QkmzN-m. (1) 

k=l 

where QkiP) = Zi{kj3) = exp(— /c/3ej) is the one-particle partition function at the tem- 
peratures kjS and Zo{i3) = 1. The microcanonical partition TnIE) can be calculated by an 
inverse Laplace-transform of (1) and is given by 

1 1 rc+ioo 

T^iE) = Y^TT- dPexpiPE)Q,{P)Zr,-M (2) 

iV ^^-j^ ZTTZ Jc-ioo 
1 N E 

= nT.I dE'r',{E')r^_,{E-E'), 



N 



where T^(E) is the inverse Laplace-transform of QkiP) and ro(-E') = 5{E). A similar, slightly 
less general, equation has recently been derived by Weiss and Wilkens ||10|| . 



Eq.(|I|) can be used to calculate all thermodynamic quantities by appropriate differen- 
tiation of InZjv. However, in any case Z^v occurs as a normalization factor and has to be 



calculated explicitly. This turns out to be a major drawback. First, the numerical effort to 
calculate Z^- grows with the square of the particle number A^. Moreover, since Z]sr{i3) grows 
exponentially with multiple precision arithmetic is required for proper calculation. We 
will present a method avoiding these difficulties. 

To ease our derivations we rewrite Zn{P) utilizing the Z-transform and define: 

Z{Z) = F{x) = E (3) 

A;=0 

Z{Q) = Gix) = £ (4) 

A;=0 ^ 

where we define QoiP) = 0. Taking advantage of the basic properties of the 2^-transform 
eq.(|l|) can be written in the form 

- x-^F{x) = F{x)G{x), (5) 



yielding [] 



F{x) = exp f f ^x-^") . (6) 



Vfc=l ^ 



Applying the inverse 2^-transform we may write Z^lf]) as 

^iv(/5) = ^/^F(x)x^-Mx, (7) 

where C := {z ^ C :\ x \= r} and r has to satisfy the condition | Z]\f{j3) \< exp{rN). 
Alternatively we may write: 



x=0 

Using eq.@ the number of particles with energy can be calculated by 



(8) 



V^{N,(3) = ^-^ In Zr,{f3) (9) 



1 ^ 



E exp{-(3kei)ZN-k{f3) 



Note that F{x) is closely related to the grand canonical partition function. 
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Some reordering yields: 

Vi{N + l,(3) = /^^fj^ exp(-/?6,)fa(iV,/3) + 1) . (10) 

Since the particle number is a conserved quantity in the canonical ensemble the direct 
calculation of the normalization factor can be omitted by using the relation 

ZN{f3) N + 1 



For Fermi systems the following recursion formula 



(11) 



V,{N + 1,P) = exp(-/3e,) (1 - r/,(iV, /?)) , (12) 



with 

ZnW) n 



(13) 



can be derived in a similar manner. In practice, only a limited number of energy levels has 
to be taken into account, since the occupation probability rapidly decreases with increasing 
energy eigenvalues. Equations (p!0[) and (|ll]) are extremely useful in practical calculations. 



The numerical effort to calculate the occupation numbers grows only linear with the num- 
ber of particles. Moreover, only a moderate arithmetic precision is required. Having the 
occupation probabilities at hand the energy expectation value is given by: 

oo 

E{N,P)=J2e^V^iN,P) . (14) 

i=0 

The calculation of the fluctuation of the occupation probabilities 6rii{N,j3) is a little bit 
more complicated and contains another recursion: 



{6v.{N + l,(3)f = l^ln(Z^.+l(/?)) (15) 



ZN+m 

+ [r/,(iV + + 1] [n,{N,(3)-ri,{N + 1,13) + 1]} 



To illustrate the usefulness of our recursion formulas we consider the ideal gas with 
parameters of liquid Helium in containers of different shapes : 
i) a cube with side-length Lx,Ly,Lz, and energy levels 

^n.,n„n. " ^^^^ ) ^'^^ 



ii) a sphere with radius a and energy levels 

2'mHeC' 



En,l = ^^ Un,l (17) 



and degeneracy an,i = 21 + 1, and 

iii) a cylinder with diameter d = 2a, height L, energy levels 

fi^ (vli m^'K^\ , . 

^-''^"'-2^^['^ + ^) ^^^^ 

with n = 1,2, . . ., I = 1,2,3, . . ., and m = . . . , —1, 0, —1, . . .. 

We denoted the zeros of the half integer Bessel functions Jn+i/2{r) by Un,i and the zeros of 
the integer Bessel function J„(r) by Vn^i- 

Fig. 1 displays the specific heats and the fluctuations of the ground state occupation 
number Srjo/N as a function of the canonical temperature for different trap geometries 
and N = 100, 1000, and 10000 He-atoms. In all cases the particle density is taken to be 
p = 0.0216A~^. With growing system size the differences between the specific heats for the 
different trap geometries almost vanish and approach the typical shape of the curve for the 
ideal Bose gas. I.e., with respect to the specific heat the boundary conditions get more and 
more unimportant with increasing volume. In contrast, the ground state fluctuations exhibit 
a completely different behaviour. The cubic box, the compact cylinder with equal diame- 
ter and height, and the sphere show almost equal ground state fluctuations for all system 
sizes, while the ground state fluctuations of the stretched box and the stretched cylinder are 
remarkably larger for temperatures below the critical temperature. This effect is not unex- 
pected because restricting the particle motion in one or two dimensions makes the system 
act like a lower dimensional system, which are known to show larger fluctuations. More- 
over, the differences between the fluctuations of the stretched traps and the compact traps 



do not decrease with increasing system size. The reason for this behaviour is found in the 
energy difference between the ground state and the first excited level, which is much larger 
for the stretched traps than for the compact traps. Since Stiq/N decreases approximately 
with A^~^/^ the infinite particle number limit is the same for all trap geometries. Under 
experimental considerations our results imply that the stability of the condensate fraction 
in anisotropic traps should be considerably smaller than in isotropic traps. In Fig. 1 we 
plotted Stiq/N to allow good comparison with previous published results Since this 

quantity goes to zero as the system size increases it is a bad indicator for phase transitions. 
The relative ground state fluctuation Srjo/riQ shown in Fig. 2 is much more conclusive in this 
respect. 

True phase transition only occur for infinite systems. However, it is well known from 



other systems, e.g. finite spin lattices and clusters |[TT|JT^ , that finite systems already display 
the onset of phase transitions. Instead of having a well defined critical temperature the 
transition occurs in a broader crossover region. As can been extracted from Fig. 2 the 
crossover region, which is indicated by the sharp increase of driQ/rjo, extends even for the 
cubic box with 10000 atoms over a temperature range of 0.5 K. For the stretched trap the 
crossover region is about twice as large. 

All calculations have been performed on an IBM-43P (233 MHz) workstation. For 1000 
particles a run with 200 temperature points took about six minutes, for 10000 particles one 
hour. A calculation of the same quantities with the recursion formula given in equation 
(HI) takes for 1000 particles about six hours (and would take at least 600 hours for 10000 
particles). Moreover, due to the numerical instabilities connected with ([l|) = 1000 was the 
largest particle number we achieved with this formula, even though we utilized a multiple 
precision package. 

We expect the recursion formulas to be quite useful for calculating the properties of dilute 
atomic Bose gases in magnetic traps with different geometry. For this propose we provide 



6 



an easy to use Java program^, which requires as input only the energy level distribution 
with appropriate degeneracy and calculates the basic properties of finite Bose systems with 
particle numbers up to 10000. A slightly faster Fortran code is available upon Email request^. 



^The code is available within the World- Wide- Web: 
http:\\www.physik.uni-oldenburg.de\ ~borrmann\BEC 

^Corresponding author: horrmann@uni-oldenburg.de 
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FIGURES 

FIG. 1. (a-c) Specific heats C^/N and (d-f) ground state fiuctuations as a function of the 
canonical temperature for systems of = 100, = 1000, and N = 10000 particles. The solid 
lines represent the results for a spherical trap, the dashed lines for a cube, the circles for a cylinder 
with a diameter to height ratio of d/L = 1, the long dashed lines for a box with side-lengths 
Lz = 4:Ly = 4:Lx, and the squares for a cylinder with d/L = 1/4. In all cases the particle density 
is taken to be p = 0.02161-^. 

FIG. 2. Relative ground state fluctuations 5r]o/r]o as a function of temperature for the cubic 
box and the stretched box and N = 1000, and 10000 particles. 
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